# Initialization ----------------------------------------------------------

library(ggplot2)
library(ggpubr)
library(data.table)
library(readstata13)
library(readxl)


# Set working directory to path of source location
setwd(getSrcDirectory(function(){}))

# Paths to the datasets and the MLE code
paths <- list(CIV = "data/civ_time_series.xlsx",
              MLE_code = "mle_cir.R")

# Read MLE code
source(paths$MLE_code)

# Read data ---------------------------------------------------------------

# Herskovi et al. (2016 JFE) CIV measure
DT_CIV <- setDT(read_xlsx(paths$CIV))[year(date) > 1945 & year(date) < 2019, 
                                      c("date", "civ")][order(date)]


# Estimation --------------------------------------------------------------

# This function is wrapper around mle_cir
# It takes as an input a dataset (DT), a variable name (varn), a time-step (dt),
# and a standardization (to know if the data needs to be divided by the mean)

est_cir <- function(DT, varn, dt, standardize = FALSE) {
  
  # Remove missing observations
  X <- DT[!is.na(get(varn))][[varn]]
  n <- length(X)
  
  # Current value is X_t and previous value is X_s
  X_t <- X[2:n]
  X_s <- X[1:(n-1)]
  
  if (standardize == TRUE) {
    X_t <- X_t / mean(X)
    X_s <- X_s / mean(X)
  } 
  
  # I use the notations 
  # dX = kappa * (alpha - X)dt + sigma * sqrt(X)dW
  mle_results <- mle_cir(X_t, X_s, Delta = dt)
  
  # Check Feller condition for stationarity
  theta <- mle_results$results
  kappa <- theta[["kappa"]]
  alpha <- theta[["alpha"]]
  sigma <- theta[["sigma"]]
  
  if (2 * kappa * alpha - (sigma^2) < 0) {
    warning("Feller condition is violated. Process is not positive a.s.")
  }
  
  return(results = mle_results)
}

# Do MLE estimation
CIR_CIV <- est_cir(DT_CIV, "civ", 1/12, standardize = FALSE)

# The object CIR_[Data] contains 4 outputs. For example, in the case of CIV
# CIR_CIV$results -- contains the coefficients
# CIR_CIV$cov -- contains the s.d. of each parameter according to the MLE
# CIR_CIV$log_likelihood -- contains the negative log likelihood at the minimum
# CIR_CIV$plot -- contains the diagnostic plots

CIR_CIV$plot
print("-----parameter estimates-----")
print(paste("stoch steady state:",round(CIR_CIV$results[["alpha"]],2)))
print(paste("mean reversion:    ",round(CIR_CIV$results[["kappa"]],2)))
print(paste("volatility:        ",round(CIR_CIV$results[["sigma"]],2)))
print("-----------------------------")




